D TV scattering and D meson resonances from lattice QCD 
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The masses and widths of the broad scalar Dp (2400) and the axial Di(2430) charmed-light res- 
onances are extracted by simulating the corresponding D-7T and D*n scattering on the lattice. The 
resonance parameters are obtained using a Breit-Wigner fit to the elastic phase shifts. The resulting 
D2(2400) mass is 351 ± 21 MeV above the spin-average \{mo + 3mn*), in agreement with the ex- 
perimental value of 347±29 MeV above. The resulting Dq -> Dir coupling g lat = 2.55±0.21 GeV is 
close to the experimental value g exp < 1.92±0.14 GeV, where g parametrizes the width V = g 2 p* / s. 
The resonance parameters for the broad Di(2430) are also found close to the experimental val- 
ues; these are obtained by appealing to the heavy quark limit, where the neighboring resonance 
Di(2420) is narrow. The calculated / = 1/2 scattering lengths are ao = 0.81 ±0.14 fm for Dtt and 
ao = 0.81 ± 0.17 fm for D*ir scattering. The simulation of the scattering in these channels incorpo- 
rates quark-antiquark as well as multi-hadron interpolators, and the distillation method is used for 
contractions. In addition, the ground and several excited charm-light and charmonium states with 
various J p are calculated using standard quark-antiquark interpolators. Our simulations are done 
in lattice QCD with two-dynamical light quarks at a mass corresponding to ~ 266 MeV. 

PACS numbers: 11. 15. Ha, 12.38. Gc 
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I. INTRODUCTION 

In the spectrum of D mesons, the IS and IP states of 
uc are well established experimentally [IHll • On the other 
hand, the knowledge of the higher radial and orbital ex- 
citations is poor, with the only experimental information 
based on a BaBar study in 2010 which found sev- 
eral new resonances whose quantum numbers are mostly 
unknown and the states are unconfirmed by any other 
experiment [l|. 

The observed pattern of masses and widths for the six 
lightest D mesons can be understood qualitatively using 
a uc valence structure and by appealing to the m c — > 00 
limit ■ The masses and widths are then independent of 
heavy quark spin s c , and the total angular momentum of 
the light quark j q = s q + L is a good quantum numbei0, 
while the total angular momentum of the state is J p = 



j p ± s c = j p ± i. The D and D* belong to the 15 
heavy quark doublet with j p = | and J p = (0 — , 1~). 
The observed narrow IP states £>i(2420) and Z?^ 2460 ) 
are candidates for j p = | + with J p = (l + ,2 + ), since 
only D-wave decays are allowed in the heavy-quark limit, 
making them naturally narrow 0- On the other hand, 
the j p = i states with J p = (0 + ,l + ) decay via S- 
wave decays in this limit Q and are related to the broad 
resonances D*(2400) and £>i(2430). 



* lmoh ler@trium f.cal 
f sasa.prclovsek@ijs.si 
^ rwww@triumf.ca 

1 We use small j when referring to heavy quark limit, while J is 
the total spin in the heavy quark limit or away from it. 



To get a quantitative description of D mesons lat- 
tice QCD can be used. In recent studies [8l-fic| masses 
for D (and D s ) mesons were calculated in lattice QCD 
with dynamical quarks close to physical values. In these 
calculations the mesons were interpolated by standard 
quark-antiquark operators. This may be problematic 
for the broad states £>5(2400) and Di(2430) which were 
not very well described by the simulation in Q. As 
well, in the charm-strange sector the scalar £>* (2317) 
and the axial D s i(2460) were discovered below DK and 
D* K thresholds, which is significantly lower than antic- 
ipated. The closeness of the masses for the scalar states 
M D » = 2318 ±29 MeV and M DsQ = 2317.8 ±0.6 MeV Q 
is not natural within a picture where the mass is domi- 
nated by the valence quark-antiquark content. This has 
led to suggestions that nearby thresholds may play an 
important role for the D* (2317) state (see, for example, 
[ill]). For these reasons a natural next step for lattice 
calculations of charm mesons is the explicit inclusion of 
multi-hadron operators and the treatment of states with 
open hadronic decay channels as resonances. In this work 
we begin with a study of D mesons where, in a lattice 
simulation, the relevant S-wave Dtt and D*tt decay chan- 
nels of scalar and axial mesons are open over a large 
range of heavier than physical pion mass. In contrast, 
a corresponding study in the D s sector would likely re- 
quire a simulation tuned to near physical quark masses 
to achieve proximity of D* and D s \ with the DK and 
D*K thresholds 0] . 

Calculations are carried out using a lattice simulation 
with dynamical u, d quarks. Correlations functions are 
constructed with quark-antiquark interpolators and, for 
the first time, also with Dtt or D*ir interpolators to study 
the relevant scattering channels. Our aim is to describe 
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the two observed broad states £>2(2400) and L>i(2430) 
as resonances, so we simulate Dtt and D*tt scattering 
and extract the corresponding phase shifts for the first 
time. A Breit-Wigner fit then allows us to extract the 
Z)q(2400) resonance mass and width, and analogously for 
the £>i(2430) within a phenomenologically motivated ap- 
proximation described in Section ITVl The remaining four 
members of the 15 and IP multiplets are stable or narrow 
in experiment. They are expected to be stable or very 
narrow also in our lattice simulation where Di(2420) and 
D 2 (2460) lie below D-wave decay threshold. We equate 
the masses of these states directly to the quark-antiquark 
energy levels on the lattice, as in all lattice simulations 
up to now. In addition, the masses of the ground and 
excited states in channels with J p = 0^,1^,2^ are ex- 
tracted. Some of these correspond to the still poorly 
known orbital and radial D meson excitations. 

The Dtt scattering in the 1 = 1/2 channel has been ad- 
dressed on the lattice only indirectly by simulating the 
scalar semileptonic D — > ir form factor fo [12] ■ Various 
scattering channels in the charm sector were simulated 
in [l3j, but the attractive 1 = 1/2 scattering Dtt or D*tt 
has not been directly simulated yet. While the scatter- 
ing lengths can not be measured, their calculation is of 
theoretical interest and we calculate the Dtt and D*tt 
scattering lengths on the lattice which can be compared 
to other types of calculations [HI, E2, IH4l7j . 

The present study of charm-light spectroscopy requires 
good control over heavy-quark discretization effects, as 
for example provided by the Fermilab method [l8| . In Q 
the spectrum of low-lying charmonium states was used to 
validate the approach. Motivated by these results on the 
low-lying charmonium spectrum a large number of non- 
exotic charmonium states up to spin 3 are studied in the 
present work. 

The present paper is organized as follows. Section [TT1 
outlines the calculational setup. Details about the gauge 
configurations, the calculation of quark propagators and 
the determination of the charm quark hopping parame- 
ter k c are discussed. Section Hill presents results for the 
spectrum of low- lying charmonium states. Encouraged 
by these results, we simulate Dtt and D*tt scattering in 
Section IIVI and extract information on scalar and axial 
resonances. For completeness some results with regular 
quark-antiquark (qq) interpolators in other J p channels 
are presented. Section [V] contains a summary and dis- 
cussion. Tables of lattice interpolating fields as well as 
details about fits and fit results are included in the ap- 
pendix. 

II. CALCULATIONAL SETUP 

Gauge field configurations were generated with n f =2 
flavors of tree level improved Wilson-Clover fermions 
fl9l . The gauge links in the action have been smeared 
using normalized hypercubic (nHYP) smearing [2lj with 
parameters (01,02,03) = (0.75,0.6,0.3). In these simu- 



Nf, x N T ki /3 a[fm] L[fm] 


#configs 


TYijx [MeV] 


16 3 x 32 0.1283 7.1 0.1239(13) 1.98 


280/279 


266(3) (3) 



TABLE I. Details of Nf = 2 gauge configurations: Nl and 
Nt denote the number of lattice points in spatial and time di- 
rections. The first error on m n is statistical while the second 
error is from the determination of the lattice scale. Observ- 
ables are based on 279 or 280 configurations. For details see 

lations the gauge fields have been generated with periodic 
boundary conditions and the fcrmion fields obey periodic 
boundary conditions in space and anti-periodic boundary 
conditions in time. The same configurations were used 
previously in a coupled channel analysis of the p meson 
|22j and in a study of K 7r-scattering [23[ • Table U lists 
some further details about the gauge configurations. For 
the determination of the lattice spacing a and the strange 
quark hopping parameter k s please refer to [22| and [23| , 
respectively. 

To calculate the quark propagation the df l_s ap_g cr 
inverter from Luscher's DD-HMC package [U, |25| | is 
used for the light and strange quarks and the same in- 
verter without low mode deflation is used for the charm 
quarks. Our final propagators are build from combi- 
nations of quark propagators with periodic and anti- 
periodic boundary conditions in time [26l . |27| . For more 
details on these so-called "P+A" propagators see [22! ]. 

A. Distillation using Laplacian Heaviside smearing 

For an efficient calculation of the quark propagation 
and flexibility in constructing correlation functions we 
use the distillation method, first proposed by Peardon et 
al. in [281 ] . In this method smeared quark sources and 
sinks are constructed using a number of low modes of 
the 3D lattice Laplacian V 2 . For an N x N matrix A 
with eigenvalues and eigenvectors one has the 
spectral decomposition 

N 

/(A) = £/(AW) (1) 
fc=i 

As in [ID, [H, [28| , the Laplacian-Heaviside (LapH) smear- 
ing with /(V 2 ) = 0(er 2 + V 2 ) is employed , so the 
smeared quark fields q s are 

N AT„ 
q s EE e^ 2 + A<*> q = £ «(*)«(*)t q , 

fc=l k=l 

(2) 

where N v depends on the target smearing a s . For this 
study we choose N v = 96 or N v = 64, depending on the 
lattice interpolating fields listed in the following Sections. 
The low mode eigenvectors and eigenvalues are calculated 
using the PRIMME package [29[. 
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B. Tuning the charm quark mass 

For the charm quarks the Fermilab method [HI, H3| 
is applied. An approach similar to the method used by 
the Fermilab Lattice and MILC collaborations [3l], [32| 
is used and we have previously employed this method to 
study charmonium and heavy- light mesons in Q. As a 
slightly modified version is used, the updated procedure 
is briefly outlined. 

In the simplest variant of the Fermilab approach 
[Hldli a single parameter, the charm quark hopping pa- 
rameter k c is determined non-perturbatively. To achieve 
this, the spin-averaged kinetic mass is measured for ei- 
ther charmonium or for heavy-light mesons and its value 
is tuned to thephysical value as determined from ex- 
periment. In |8[ the spin- average of the IS states in 
the spectrum of D s mesons was used for this purpose. 
As we here use gauge configurations with only 2 dy- 
namical quark flavors we instead opt for the IS char- 
monium states and tune the spin-averaged kinetic mass 
(M Vo +3Mj/y)/ 4 to its physical value. As further param- 
eters, the Fermilab action contains the clover coefficients 
ce and cb and, incorporating tadpole improvement we 
choose ce = cb = ciw = 4r where uq is the average link. 
To determine the average link we calculate the Landau 
link on unsmeared gauge configurations. This leads to 
the numerical value ciw = 1.75218 for the ensemble used 
in this study. 

With the description outlined above, the remaining 
task consists of determining the kinetic mass M 2 by em- 
ploying the general form of the lattice dispersion relation 
from |32j 

yp ' 2M 2 6 8M? ' w 

i * 

where p = for a given spatial extent L. Even ne- 
glecting higher orders and taking only the terms explic- 
itly listed, this form contains too many parameters. We 
therefore determine Mi using two simplified methods: 

(1) neglect the term with coefficient W4 and fit Mi, 
M 2 and M 4 . 

(2) fit E 2 (p) and simplify the (p 2 ) 2 term arising from 
the mismatch of Mi, M2 and M4 by setting Mi = 
M4 for charmonium and M 2 = M4 for heavy-light 
mesons. 

Note that (2) differs slightly from the method previ- 
ously used [8j. This change is motivated by the results 
from method (1) where we obtain Mi ss M4 for char- 
monium and M2 « M4 for heavy-light mesons. The 
modified method (2) therefore tests if fits qualitatively 
change when including the term — - ^ 4 ^ pf breaking 
rotational symmetry. 

For the tuning of k c we used correlation functions with 
sources on 2-4 time slices on the full ensemble. A cross 



check with heavy-light and heavy-strange states used 
data on 16 time slices. Tables HU IIIII and IIVI list the 
results for our final choice k c — 0.123. 





Method (1) 


Method (2) 


Mi 


1.52499(42) 


1.52484(42) 


M 2 


1.9581(59) 




M 4 


1.5063(216) 




M 2 
M t 


1.2840(38) 


1.2745(41) 


M 2 [GeV] 


3.1186(94)(327) 


3.0951(102)(325) 


Exp [GeV] 


3.06776(30) 



TABLE II. Fit parameters obtained for charmonium with 
both tuning methods. The values in the last two rows are 
in GeV, while all other values are in lattice units. The first 
error on the kinetic mass M2 is statistical while the second 
error is from the scale setting. The results for M4 are not 
used in our setup but displayed to demonstrate the observed 
relation Mi « M4. The last row contains the experimental 
value from [33| . 





Method (1) 


Method (2) 


Mi 


1.09704(57) 


1.9716(58) 


M 2 


1.2917(61) 




Ma 


1.2397(144) 




M 2 
Mi 


1.1774(56) 


1.1771(65) 


M 2 [GeV] 


2.0572(98)(216) 


2.0568(116)(216) 


Exp [GeV] 


2.07634(38) 



TABLE III. Same as Table [TjJ but for charm-strange (D a ) 
mesons. 

In all three cases methods (1) and (2) are in reason- 
able agreement. As a subset of our charmonium data 
has been used for tuning k c , it is no surprise that the 
kinetic mass of the charmonium spin-averaged ground 
state agrees reasonably well with the PDG value. The 
spin-averaged D s mass also compares favorably to the 
experimental value. While the charmonium result sug- 
gests that our charm quark mass has been tuned just 
a tiny bit to heavy, the result for the heavy-light system 
comes out somewhat lower than expected for our unphys- 
ical light-quark mass. The heavy quark discretization ef- 
fects for heavy-heavy and heavy-light systems differ and, 
as we are dealing only with the simplest (lowest order) 
version of the Fermilab action a perfect agreement is not 
expected, especially on a rather course lattice. We how- 
ever conclude that our quark masses are tuned reasonably 
well for the current purpose. This will be tested further 
in section IIIII where the low- lying charmonium spectrum 
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Mpthnd (TI 


Mpthnd (2) 

1VJ.C L11U H \ _ J 


Mi 


1.04226(111) 


1.04206(113) 


M 2 


1.2242(161) 




M 4 


1.2550(461) 




M 2 
Afl 


1.1745(154) 


1.1573(168) 


M 2 [(?eV] 


1.9497(257) (205) 


1.9207(288) (202) 


Exp [GeV] 


1.97140(13) 



TABLE IV. Same as TableUQbut for charm-light (D) mesons. 
The results for M4 are not used in our setup but displayed 
to demonstrate the observed relation M2 ~ M4. Notice that 
the value for M2 in physical units is based on a heavier than 
physical light-quark mass. 



is calculated. 

C. Variational method and correlator basis 

To extract the low-lying spectrum we calculate a ma- 
trix of correlators at every source and every sink time 
slice ti and tf 

C»(t = t f - U) = J2(0\Oi(tf)O}{U)\0) (4) 

ti 

= J2^ tEn (m\n)(n\O]\0), 



using suitable lattice interpolating fields with definite 
quantum number^ J PC (for charmonium) or J p (for 
heavy-light states). To extract the low- lying spectrum 
the generalized eigenvalue problem is solved for each time 
slice 



C{t)^ = \ {k \t)C(t )J {k) , 
A (fc) (i)«e'* (l + 0(e"* A ^ 



(5) 



)) 



At large time separation only a single state contributes 
to each eigenvalue. This procedure is known as the vari- 
ational method j33 - l36j . The employed interpolators are 
listed in (|21I23|I and the Appendix [5] 



III. CHARMONIUM RESULTS 

Recent lattice simulations of excited charmonium 
states were presented in [H, fic l W H, |38| . The mixing of 
cc and DD was explored in [37|. higher spin and ex- 
otic J PC states with carefully determined continuum spin 
were presented in j38| . 



In our previous study [8| the low-lying charmonium 
states have been determined as a benchmark for our 
heavy-quark treatment. The resulting low- lying spec- 
trum was in good qualitative agreement with experiment. 
For our current study, the distillation method based on 
lowest eigenmodes of the lattice Laplacian enables us to 
have considerably more freedom with regard to the lat- 
tice interpolating fields used. We exploit this opportunity 
and use an enlarged basis which is tabulated in Appendix 
IA II We consider all non-exotic channels up to spin 3 and 
also aim at extracting excited energy levels in channels 
where these are expected above multi-particle thresholds. 




Al + Tf" A1 ++ T1 ++ T2 ++ T2 + T2" A2" A 2 + A2 ++ 

FIG. 1. Energy differences AE = E - j(M Vc + 3M. H ,) for 
charmonium states on the lattice and in experiment; refer- 
ence spin-averaged mass is j(M Vo + 3Mj^,) w 3068 MeV in 
experiment. Levels are listed according to lattice irreducible 
representation. The results from simulation are displayed 
as red stars and displaced slightly to the left, while our new 
results displayed as blue crosses are displaced slightly to the 
right. The statistic and scale setting uncertainties have been 
combined in quadrature. Experimentally observed states are 
plotted as black bars or (where there is a substantial uncer- 
tainty in mass determination) filled boxes. The level corre- 
sponding to the well established X(3872) has been plotted 
in both the Tl ++ and T2~ + irreps, reflecting its uncertain 
quantum numbers [39l , |40T | . 



As usual J is the spin, P is parity and C charge conjugation. 



Figure [T] collects results for all lattice irreps and non- 
exotic quantum numbers. For all states the difference 
with respect to the spin-averaged ground state Mjg = 
(M Vc + 3M//$)/4 is plotted. The use of a larger ba- 
sis enables us to extract at least the ground state in all 
channels investigated and in many cases also one or more 
excited state energy levels can be extracted. Details of 
the fitting methodology can be found in the appendix. In 
addition to our new data (blue crosses) we also display 
the results from [1] (red stars). 

Figure [2] shows the same energy levels as before, as- 
signed to continuum states. For the assignment, degen- 
eracies across irreps and interpolator overlaps have been 
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usecjf], including data from the E irreducible representa- 
tion which is not shown in Figure [1] 
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FIG. 2. Energy differences AE = E - \(M Vc + 3Mj^) for 
charmonium states on our lattice and in experiment; refer- 
ence spin-averaged mass is \{M Va + 3Mj^,) w 3068 MeV in 
experiment. In addition to the low-lying IS, IP and 2S states, 
the full set of ID and 2P states is also seen. The black lines 
or boxes denote known experimental states with uncertainty 
on their masses. The level corresponding to X(3872) is plot- 
ted for both possible quantum numbers 1 ++ and 2~ + . The 
magenta lines on the right denote relevant lattice and contin- 
uum thresholds. From low to high energy, the physical DD, 
D*D and D* D* thresholds are plotted as dotted lines. The 
corresponding lattice thresholds are plotted as dashed lines 
in a similar order: DD in S-wave, D*D in S-wave, DD in 
P-wave and D* D* in S-wave. The corresponding D s meson 
thresholds are omitted as calculations were performed with a 
2-flavor sea. 

The dotted lines in Figure [2] denote the physical DD, 
D*D and D*D* thresholds. Our results for states below 
all thresholds agree well with the experimental states, 
which are commonly interpreted as the IS, IP and 25 
multiplets (from low to high mass). Above the physi- 
cal DD threshold we observe another band of states of 
quantum numbers 1 , 2 , 3 and 2 ^ which is nat- 
urally interpreted as the ID states. Notice that we seem 
to observe the D-wave 3 state also in the Tl~~ irrep. 
The J PC = 1" D-wave state corresponds to the experi- 
mental 'I' (3770) which has an appreciable decay into two 
D mesons. Note that this decay proceeds in P-wave and 
the corresponding lattice threshold is far away from the 
physical DD threshold. As a result we do not expect to 
reproduce the correct \&(3770) mass or the correct split- 
ting between the ^(2S) and 'I' (3770). Similar remarks 
can be made for all other states above threshold, so we 
restrict ourselves to some qualitative observations. Of 
particular interest is the observation of a further band of 



states split from the spin averaged 15* states by around 
850 MeV in our simulation. These states have the pat- 
tern expected for the 2P states. Notice that the Xc2{2P) 
has been identified in experiment and the corresponding 
energy level is included in our plot. Furthermore, a recent 
study by BaBar [43| suggests that the A(3915) is likely 
to have quantum numbers ++ and could be interpreted 
as the Xco(2P)- The identification of this state with the 
X c o(2P) is not appealing on theoretical grounds Q and 
there arc indications for a broad Xco{2P) [44| at a lower 
mass, which would be more compatible with our dat£@. 
Evidence for the ^2 with quantum numbers J p = 2~~ 
has furthermore been found by the Belle collaboration 
@ at a mass of 3823.5(2.8) MeV. With regard to the 
X(3872), whose quantum numbers are not yet settled 
[39l Eoj and could be either 1 ++ or 2 we find that 
one of the states we interpret as a 2P state is close to the 
mass of the X(3872). Just like in experiment, this state 
is compatible in mass with the D*D threshold. It should 
be stressed that the current results are within uncertain- 
ties compatible with both possible quantum numbers and 
that one can not draw any strong conclusion about the 
nature of the X(3872) from this study. Nevertheless it is 
interesting that a state is observed in close vicinity to the 
threshold in the 1 ++ channel, while the ground state in 
the 2 h channel comes out lighter. For a recent lattice 
study investigating this issue see [Irjj . For a particularly 
insightful discussion regarding the possible nature of the 
X(3872) see [13]. 

In addition we observe a number of further states for 
which likely assignments are shown in the Figure. In 
particular we find two spin 3 F-wave states and another 
set of excited S- waves. 

To disentangle spin-dependent from spin-independent 
contributions we further define spin-averaged masses 

M ss =±(M v , a + 3M J/v ), 

M w = i (M Xc0 + MI Xcl + 5M Xc2 ) , (6) 

M2P = l(M x , o +3M x , i+ 5M x , 2 ) . 

The results are listed in Table [V] In addition we take a 
look at the hyperfine splittings between spin-singlet and 
spin-triplet states 

M„ 3i - M n i L (7) 
and at the P-wave spin-orbit and tensor splittings 

Mspm-orut = g(5M Xo2 - 3M Xol - 2M Xc0 ) , (8) 

MTensor = £ (3M Xcl - M Xc2 - 2M Xc0 ) . 



3 For a more elaborate way to identify the continuum spin please 
refer to SE El 



4 For the splitting between the Xc2(2P) and the Xco(2P) we obtain 
63 ± 33 MeV. 
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Depending on the heavy quark treatment, lattice dis- 
cretization effects in these quantities can be substantial. 
Their determination is a challenge for lattice QCD. Val- 
ues are extracted for the 15, 25 and IP hyperfine split- 
tings and for the IP and 2P spin-orbit and tensor split- 
tings and are presented in Table IV1 The experiment val- 
ues in the table are the corresponding PDG values [33j . 
In the case of the 15 rj c state which enters the hyperfme 
splitting, the PDG average has a poor confidence level 
and newer results from BESIII [48| and Belle [4£| sug- 
gest that the hyperfme splitting is substantially lower. 
For recent lattice results on the 15 hyperfme splitting 
including a continuum extrapolation see j50l - l52j . 

For the charmonium hyperfme splitting we also deter- 
mine the uncertainty associated with the kappa tuning 
procedure outlined in Section III Bl First we average over 
the results from method (1) and method (2) for our tun- 
ing runs at k c = 0.123 and k c = 0.124 which are close to 
the value corresponding to physical Mi. We determine 
the resulting values for Mi and the hyperfme splitting for 
both values of k c . Due to the enlarged statistics our final 
charmonium data differs slightly from the tuning run at 
the same k c and we also determine the kinetic mass Mi 
for our final data. We then use interpolations of the tun- 
ing data to determine the kappa tuning uncertainty for 
our final data. For this purpose the uncertainty from our 
choice of fitting model is estimated by the difference be- 
tween the kinetic masses obtained from method (1) and 
method (2). For the total error, the stochastic error from 
the Monte Carlo estimation, the scale setting uncertainty 
and the uncertainty from the fitting model are added in 
quadrature. As our final run turns out to have slightly 
missed the physical Mi we obtain an asymmetric error 
of ±o'q from the uncertainty in the charm quark mass. 
Similar tuning errors of about 2% are expected for all 
spin-dependent quantities. 

In section Hi Bl we also took a look at the kinetic masses 
of D and D s mesons. In addition to charmonium mass 
splittings the values for the combinations 2M-g—M^ and 
2M-jj — Mjg are also provided in Table fVl In these combi- 
nations the leading heavy quark contribution drops out. 
Again the proximity of our results to the experimental 
values is encouraging. 



IV. D MESONS RESONANCES INCLUDING 
MESON-MESON INTERPOLATING FIELDS 

A. Energy levels and the Liischer method 

This section provides a short but general introduction 
for the extraction of resonance parameters. A small mod- 
ification specific to our heavy-quark setup is discussed in 
the following section. 

Assuming a localized interaction, the energy levels of 
a two-hadron system in a finite box are related to the 
scattering phase shift in the elastic region [551-58] . One 
first needs to determine energy levels E of the two-hadron 



Mass difference 


This paper [MeV] 


Experiment [MeV] 


Xco — IS 


354.7 ±4.4 ±3.7 


347.0 ± 0.4 


Xci - is 


425.7 ±3.9 ±4.7 


442.9 ±0.3 


Xc2 - IS 


468.7 ±4.7 ±4.9 


488.4 ±0.3 


h c -IS 


438.0 ± 4.9 ± 4.6 


457.7 ±0.4 


TP -IS 


441.7 ±4.0 ±4.6 


457.5 ±0.3 


v'c-Ts 


548.9 ±4.9 ±5.8 


569.2 ±4.0 


J/*' - IS 


606.8 ±4.9 ±6.4 


618.3 ±0.3 


2S -IS 


592.3 ±4.9 ±6.2 


606.1 ± 1.0 


IS hyperfme 


107.9 ± 0.3 ±l.l±ao 


116.6 ± 1.2 


IP spin-orbit 


39.7 ±2.1 ±0.4 


46.6 ± 0.1 


IP tensor 


11.02 ±0.87 ±0.12 


16.25 ±0.07 


IP hyperfine 


3.7 ±2.7 


-0.10 ±0.22 


2S hyperfine 


57.9 ±2.0 


49 ±4 


2P spin-orbit 


24.6 ±15.7 ±0.3 


- 


2P tensor 


2.2 ±4.3 




2P-TS 


836.4 ± 30.5 ± 8.8 




271% - M- c 


890.9 ±3.3 ±9.3 


882.4 ±0.3 


2M-JJJ - M- c 


1065.5 ±1.4 ±11.2 


1084.8 ± 0.6 



TABLE V. Mass differences in the charmonium spectrum in 
MeV compared to experimental values (calculated from [33f ] : 
the value for the IP hyperfine splitting is from [53j]). Bars de- 
note spin-averaged values. For the results of this paper, the 
first error denotes the statistical uncertainty and the second 
error denotes the uncertainty from setting the lattice scalfl 
In addition there is a non-negligible error from the uncertainty 
in the determination of n c for all spin-dependent quantities. 
For the IS hyperfine-splitting the corresponding error is esti- 
mated and given as the third (asymmetric) error. It is stressed 
that the gauge ensembles at our disposal do not allow for a 
continuum and infinite volume extrapolation. Consequently 
qualitative but not quantitative agreement is expected. In the 

last line we also provide the splitting 2A% A%j which can 

be directly compared to the results quoted by the Fermilab 
lattice and MILC collaborations [3lJ and also to the value of 
2M Ds - M Vc quoted by HPQCD in [H]. For the determina- 
tion of the strange quark mass on our lattices please refer to 

E3. 

a For spin-dependent quantities the indirect contribution of the 
scale setting uncertainty to the kappa tuning uncertainty is 
sizable. Our scale setting error only accounts for the direct 
uncertainty associated with the setting of the lattice scale. 



system on the lattice. We choose the total momentum 
of the system to be zero in this simulation, so the lattice 
frame coincides with in the center of momentum (CM) 
frame, and both hadrons have momentum p*. In this 
case we avoid the complications that arise for the extrac- 
tion of the S-wave from simulations with non-zero total 
momentum caused by ^ m w [o9l.l60|. In the exterior 
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region, where the interaction is negligible, 



E z = s 



(9) 



yr 2 + m 2 



Hi 



+ 



j^ 2 + m 2 



and the discrete values of p* are extracted from E via 

p* 2 = — (s - (m Hl + mH 2 ) 2 ) (s - (m Hl - m H2 ) 2 ) , 

(10) 

while the corresponding dimensionlcss momentum q is 
defined as 



q= 2^ P 



(11) 



Determining the value of the momentum q from these 
relations, one obtains the relevant S-wave scattering 
phase shift from the Liischer formula [57| 



tan<5(<7) = 



7T2 q 



Zoo(l;g 2 ) 



(12) 



where Zqq(1; q 2 ) is a generalized zeta function. This for- 
mula neglects terms exponentially suppressed with the 
lattice volume, and we note that these terms might not 
be completely negligible for our volume. We are setting 
up for simulations at a larger volume and results will 
indicate whether this might affect our present results. 

The elastic scattering phase 8 is related to the scatter- 
ing amplitude TJ by 



Ti = sin<5 ; e lSl = 



2i 



A variable pi(s), defined as 

(p*) 2 ^ 1 cot S(p*) 

Pl{s) = 75 ' 

depends on the scattering length ai near threshold 

V~spi(s) = -+0(p* 2 ). 
ai 



(13) 



(14) 



(15) 



For the case of an elastic channel dominated by a single 
resonance one can also assume a relativistic Breit-Wigner 
shape and obtain 



y/~sT r (s) cot S(q) = s r 
T r (s) = g~ 



(16) 



where the width of the resonance T r has been 
parametrized by a phase space factor and a coupling con- 
stant g; in our case g will either be gD'Dir or gDxD*-*- We 
will extract the mass m r and the coupling g from the re- 
sulting pi(s) via 



Pi(s) = -4(s r - s) 

r 

obtained by combining Eqs. (|16[) and ([1 



(17) 



1.35- 
1.3 

1.25 
1.2 

1.15 
1.1 

1.05 
1 

0.95 




3 4 
q 2 =(Lp*/27t) 2 



FIG. 3. Measured energies of D or D* as a function of 
q 2 = (Lp* /2tt) 2 and a fit using the dispersion relation © 
with method (1): Mi » 0.980, M 2 » 1.107, M 4 ~ 1.107 are 
extracted for D, while Mi « 1.063, M 2 w 1.267, M 4 ~ 1.325 
are extracted for D* . 



B. Dispersion relations 

As already mentioned we have to modify the above 
procedure slightly to account for our unphysical disper- 
sion relation. Equation Q uses a relativistic dispersion 
relation for both hadrons. In our case one hadron is a 
pion, where the relativistic continuum dispersion relation 
holds well [23| in the simulation. The other hadron is a 
heavy-light D or D* meson with the dispersion relation 
([3]) displayed in Fig. [3l where the values for Mi , M-i and 
M4 are determined separately for D and D* with method 
(1) and provided in the caption. Therefore, we extract 
the momentum p* from the energy E of the D^*'ir system 
via 



Mi 



2M 2 



8M| 



(18) 



For convenience the values of separate energies for tt, 
D and D* at momenta q = and q = 1 are listed in 
Table I VII Their sums provide the reference energies of 
the lowest two non-interacting scattering states and are 
plotted as dashed lines in Figs. |3]and[7] 



C. Results 

For the D mesons a basis consisting of quark-antiquark 
and meson-meson interpolating fields is used. For the 
qq part the interpolating fields are tabulated in Table 
Kill For mesons made from quarks with different masses, 
charge conjugation (or more generally G-parity) is not a 
good quantum number. In case of the J p = 1 + D\ meson 
we therefore consider also mixing between interpolating 
fields corresponding to different charge conjugations in 
the mass-degenerate case. This has already been found 
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meson 


Ea (p* = 0) 


E dr a (p* = ^) 


71 


0.1673(16) 


0.4268(65) 


D 


0.9801(10) 


1.0476(10) 


D* 


1.0629(13) 


1.1225(14) 



TABLE VI. Energies of tt, D and D* for momenta p* = 
, which are relevant for non-interacting D'*^ (0)-7r(0) and 
(l)7r(— 1). The pion energy at p* = 3j£ is based on the 
continuum dispersion relation; the energies of D and D* are 
obtained with the dispersion relation © within method (1) 
and with Mi,2,4 from Fig. [3] 



to be important in Q and has been investigated for kaons 



Gl] 



For the I = ^ resonances in the J p = + (Dq) and 
J p = 1+ (£>i) channels, one needs the following meson- 
meson combinations 



1 1 

2' 2 



- 1 ) 

2 2/^o 



D~7T 



D°7T° 



(19) 



D*-7r+ + \ -D*"tt 



with 



£T = grd , (20) 
7r + = dTu , 

7T° = — y= (uTu — dTct) , 

and r = 75, 7^ for pseudoscalar and vector fields, respec- 
tively. In case of the Dq we use a basis of six interpolating 
fields. The first four are qq interpolators as listed in Table 
IXIII and the last two are meson-meson interpolators: 

5 (t) = ^D-(0)n+(0) + y|^°(0)7r°(0) , (21) 



D-( Pl )n + (-p t ) + ^ ( Pl h°(-P* 



Interpolator Og is constructed from nontrivial momenta 
^Pi with 



Pi 


= (1,0,0), 


Pi 


= (-1,0,0), 




P2 


= (0,1,0), 


Po 


= (0,-1,0), 


(22) 


P3 


= (0,0,1), 


P6 


= (0,0,-1). 





For the D\ a basis of ten interpolating fields is used. 
Again just two of these arc meson-meson interpolators 
constructed in a similar way 



09 (*) = 

0ioW = E 



-(0)7r+(0) + ^-5* (0)7r>), (23) 
Hi)*-(p i ) 7 r+(-p i ) + J±D*° ( Pi )7r (- Pi ) 



1.6 



^1-4 



1.2 



0.8 



H m 



"DfflMSLS. 



« a » 



12 



15 



_J 

18 



21 



FIG. 4. Lowest three energy levels in the Dg channel with 
J p = + . The fit ranges and uncertainties are indicated by 
the solid horizontal lines. In addition, non-interacting scat- 
tering levels D(p)n(—p) are depicted by dashed lines. Inter- 
polator choices and numerical values can be found in Table 

rvm 



In case of the interpolating fields with non-trivial mo- 
menta, we restrict the number of Laplacian eigenmodes 
used in the construction to N v = 64, while N v = 96 is 
used for all other interpolating fields. The contractions 
for such 7=1/2 interpolators are explicitly provided in 



23(, with 



Appendix B of the Ktt scattering simulation 
the only necessary replacement s — > c. 

The fitting methodology is the same as for charmo- 
nium and is outlined in Appendix IA 21 Table IXIIII in 
Appendix IA 31 lists our choices for the interpolator sets, 
timcslice to and fit ranges as well as the fit results and 
% 2 /d.o.f. for the fit in the channels without meson-meson 
interpolating fields; the basis used is indicated by the in- 
terpolator numbers from Tabic IXlTl The energy levels for 
these channels are collected in Figures [10] and [TT] The 
results for the Dq and D\ channels that take into account 
the meson-meson interpolators are discussed separately 
below. 



1. Dtt scattering in the J = + channel and Dq resonance 

We now consider the J p = + channel where only one 
resonance has been established in experiment [Hj], the 
£>o(2400). The first step is to extract the energy levels 
in the finite volume. Figure [4] shows the effective masses 
obtained for the three lowest states in the Dg channel for 
our final selection of interpolating fields. Results from 
the full basis agree qualitatively but are more noisy. For 
all displayed states, we can obtain stable fits which are 
reasonably insensitive to the choice of fit range, number 
of fit exponentials (one or two) and choice of to . 

To illustrate the effect of a combined basis consisting of 
both quark-antiquark and meson-meson interpolators we 
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level n 


interpolators 


to 


fit range 


type 




X /d.o.f 


ap* 


„„2 
sa 


do 


1 


1,3,5,6 


3 


4-21 


2 exp 


1.1145(25) 


3.39/13 


0.0939(34)i 


1.2420(65) 


41.2(12.2)i 


2 


1,3,5,6 


3 


4-13 


2 exp 


1.3060(52) 


4.73/6 


0.2474(51) 


1.7057(135) 


-77.1(2.8) 


3 


1,3,5,6 


3 


4-11 


2 exp 


1.495(15) 


0.35/4 


0.4093(118) 


2.236(45) 


-16.7(12.4) 



TABLE VII. Final results for the lowest three energy levels in the Dq channel with J p = + . For each state the timeslice to 
for the variational method, the fit range, fit type and x 2 /d.o.f as well as results for the energy E n , the momentum p* defined 
in (|10|l . the invariant mass squared s and the S-wave scattering phase So are provided. Interpolators O1-4 of type uc are listed 
in Table IXlTl while 05,6 of type Dn are given in (|21|l . 



level n 


interpolators 


to 


fit range 


type 


E n a 


X 2 /d.o.f 


ap* 


so 2 


So 


1 


1,4,7,8,9,10 


3 


9-19 


1 exp 


1.1978(28) 


7.30/9 


0.0938(40)i 


1.4348(68) 


40.9(14.0)i 


2 


1,4,7,8,9,10 


3 


8-16 


1 exp 


1.3222(90) 


1.72/7 


0.1810(291) 


1.748(24) 


-43.3(4.7) 


3 


1,4,7,8,9,10 


3 


8-16 


1 exp 


1.3456(71) 


5.34/7 


0.2068(80) 


1.811(19) 


-55.7(4.0) 


4 


1,4,7,8,9,10 


3 


7-11 


1 exp 


1.571(10) 


0.30/3 


0.4107(83) 


2.469(32) 


-16.1(7.0) 



TABLE VIII. Final results for the lowest four energy levels in the Di channel with J 
see Table Ell] 



1 + . For an explanation of the entries 



qq: O I Dtc: (), ,, 



just qq: () , . 



justDlt: O 
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FIG. 5. Energy levels in the J p = + channel for three differ- 
ent choices of interpolator basis. The panel on the left shows 
the full results from a basis d,3,5,6 of uc and Dtv interpo- 
lators. The mid panel shows results from just uc interpola- 
tors (01,3,4), while the right panel from just D7r interpolators 
(05,6)- Interpolators are listed in Table IXlTl and Eq. (f2l| . 
All data are for to = 3. 



compare the results for the lowest levels in the Dq chan- 
nel for different choices of interpolator basis in Figure [5] 
The results plotted in the left panel correspond to the 
data from our final choice of interpolators already shown 
in Figure 2J In the right panel we plot data from our 
2x2 basis of meson-meson interpolators. As expected 
these interpolators lead to energy levels in the vicinity 
of the non-interacting D(0)n(0) and D(1)tt(— 1) states, 
but also show sizable excited state contaminations. The 
mid panel shows results from a basis consisting of a sub- 
set of qq interpolators. The ground state for this choice 



0.2 







CO 



-0.2 
CJ '-0.4 



o 




* 

Sb -0.6 



1.2 



1.4 



1.6 



1.8 2 
s 



U 1 I 

2.2 2.4 



2.6 



FIG. 6. The quantity po(s) = cot S for the Dir scattering in 

channel J p — + as a function of s, both in units of the lattice 
spacing. For a single Breit-Wigner type resonance points are 
expected to lie on a straight line as suggested by (|16I17[) . 



turns out to be in the vicinity of the D(0)n(0) state but 
has much larger errors than the ground state from the 
full basis in the plateau region. The effective masses 
calculated from the second and third eigenvalues never 
plateau and are very noisy. This is quite contrary to the 
full basis, where the plateau for the level n = 2 is well 
determined. From this plot it is quite obvious that an 
analysis of the energy levels considering only qq interpo- 
lating fields would not lead to satisfactory results with 
our sources and statistics. 

Table [VTT] shows the results for the preferred interpola- 
tor choices that combine 1717 and meson-meson interpola- 
tors and correspond to the levels in Figure |U It provides 
the momentum p* defined in (|10[) . the invariant mass 
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squared s and the S-wave scattering phase So for the three 
lowest levels. The ground state energy in this attractive 
channel is below the non-interacting D(0)n(0) level and 
the corresponding p* and Sq are imaginary. While the 
phases for the first two levels are fairly well determined, 
our conservative estimate for the third level differs only 
fairly insignificantly from the non-interacting D(l)ir(— 1) 
level. As a consistency check, we therefore compared our 
results for the energy levels with values from the ratio 
method, used, for example, in the case of the p meson by 
the PACS-CS collaboration [HIH. Within errors, the 
extracted energy levels agree. 

The S-wave Dir scattering length ao = tan(<5)/p* (fT5|) 
is extracted from the lowest level with small p* 



a~ 1/2 = 6.56 ± 1.16a 



0.81 ± 0.14 ± 0.01 fm , 



(24) 



7=1/2 
^0 



= 17.7 ± 3.1 ± 0.2 GeV~ 



where ^ from the simulation were inserted to the re- 
duced mass hdtt- The ratio ao/pp^ is independent of 
within Weinberg's current algebra result aplpp-K = 
1/(2ttF 2 ) ps 9.4 GeV~ 2 with F n ps 0.13 GeV [ll, [T3l- 
Heavy meson ChPT combined with the lattice input from 
[13 leads to ao/^DTT ~ 9 - 11 GeV" 2 [HG3 with P h y si - 
cal m,D,TT input to the reduced mas^|. An indirect deter- 
mination based on the simulation of the D —> ir semile 
tionic transition leads to aa/pD-K = 16.4±2.3 GeV -2 
which agrees with our result within error. A calculation 
using Unitarized ChPT and taking into account coupled 
channel effects results in a$/ pd-k = 13.8 ± 0.4 GeV -2 
which is also compatible with our result. It is interesting 
to note that our a /p ~ 18 GeV -2 for Dir is very close 
to our result for D*ir ^7} and Kir [H| in I = 1/2 chan- 
nels. Indeed, current algebra predicts the same ratio for 
all three channels, albeit the current algebra result itself 
is lower. 

In Figure El we plot p (s) = ^cot<5 (0 for the D* 
channel as a function of s. For a single Breit-Wigner 
type resonance a linear relationship (|17[) should emerge. 
Unfortunately this relationship can not be tested with 
our current data, as our highest energy level is not de- 
termined well enough. However, in the elastic regime 
where the phase shifts 6q are rising monotonically, our 
data for the two lowest energy levels clearly indicates a 
resonance in between levels 1 and 2. Assuming a Breit- 
Wigner shap^l we extract the resonance mass and cou- 
pling g l £i Dyr from po(s) with a linear fit (jTT)) over three 



5 Where given we use the masses provided by the authors for cal- 
culating the reduced mass. When not provided we use the values 



from [16 

Note also that all phases have been extracted under the assump- 
tion that admixtures from higher partial waves due to the broken 
rotational symmetry on the lattice do not play a significant role. 
In the case of the frame with total momentum 0, these admix- 
tures enter only at I > 4 and should be small. 



pointsQ. We obtain 



Aat 



lat 
9d*Dtt 

ji flat 



2.55(21)(03)GcV , (25) 
350.8(20.2)(3.7)MeV , (26) 



where Mjg = \{Md + 3Mc»). The resulting mass is 
compared to the PDG value for the Dp (2400) m Figure 
1101 To compare our coupling to experiment, we can use 
the total width T D * = 247(40) MeV [H] and translate it 
into an upper bound for the coupling. The resulting value 
9d*d-k — 1-92(14) GeV is not too far from our estimate, 
although slightly smaller. 



2. D*tt scattering in the J =1 
resonances 



channel and Di 
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7. Lowest four energy levels in the Dq channel with 
= 0+. The fit ranges and uncertainties are indicated by 
the solid horizontal lines. In addition, non-interacting scat- 
tering levels D* (p)tt(— p) are depicted by dashed lines. Inter- 
polator choices and numerical values can be found in Table 

lYnn 

Unlike the previous case, there are two known reso- 
nances m the J p = 1+ channel, the D,(2420) and the 
Di(2430) [33|]. Again, we start by extracting the energy 
levels for the lowest states from our simulation. The re- 
sults are plotted in Figure [7] From the experience with 
the Dq in the previous section, we expect to extract two 
energy levels in the vicinity of the lowest scattering states 
-D*(0)7r(0) and D*(l)ir(— 1) and two additional levels re- 
lated to the two resonances in this channel. As can be 
seen in Figure [7] this is the case. 

Again one can compare the results from the final choice 
of interpolators to subsets containing only qq or only 



7 Due to the large error for po(n = 3), the results are almost 
independent if the level n = 3 is taken in the linear fit or not. 
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FIG. 8. Energy levels in the J = 1 + channel for three differ- 
ent choices of interpolator basis. The panel on the left shows 
the full results from a basis ©1,4,7,8,9, 10 of uc and D*n interpo- 
lators. The mid panel shows results from just uc interpolators 
(0i,4,7,s), while the right panel contains our results from just 
D*tv interpolators (09,10). Interpolators are listed in Table 
EH and Eq. (|23]l. All data are for to = 3. 
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FIG. 9. Same as Figure [6] for the Di channel with J p — 1 + 



meson-meson interpolators. This comparison is shown 
in Figure |S] The results from qq interpolators alone are 
shown in the mid panel. The two largest eigenvalues lead 
to effective masses which seem to display a clear plateau 
at intermediate source-sink separations. In the right 
panel the results using only meson-meson interpolators 
are shown. Here we observe clear signals in the vicinity 
of the non-interacting D*(0)7r(0) and D*(l)7r(— 1) states. 
Notice that the lowest level from just qq interpolators 
is at best marginally compatible with the ground state 
from meson-meson interpolators. Turning our attention 
to the full basis shown in the left panel, we notice that 
the ground state is compatible with the state observed 
from meson-meson interpolators alone, while the n = 2 
level is very similar to the n — 2 level with just qq in- 



terpolators. There is no level in the vicinity of the qq 
ground state (green left triangles in the mid panel) but a 
new state (blue diamonds in the left panel) emerges. The 
n = 4 state is found in the vicinity of the non-interacting 
£>*(l)7r(— 1) level. It is interesting to see that one of the 
levels observed with just qq interpolators survives with 
no significant change in energy, while the lowest changes 
quite drastically. We will return to this observation for 
our interpretation of the data below. 

The numerical results and fit parameters correspond- 
ing to the final choice of basis are tabulated in Table 
IVIIII For the three lowest levels the S-wave phase shift 
do is well determined. For the fourth level, which is in 
the vicinity of the non-interacting D*(1)tt(— 1) state and 
has large overlap with interpolator O10, the results have 
a large uncertainty, just like in the case of the third level 
in the J p = + channel. The results for po(s) = 2j= cot Sq 

n 

arc plotted as a function of s in Figure \9\ 

The S-wave D*tt scattering length (fT5|) is extracted 
from the lowest level with small p* 



1=1/2 



a Q • =6.53 ± 1.34 a 



= 0.81 ± 0.17 ± 0.01 fm 



(27) 



a 



1=1/2 




P-Dtt 



17.6 ± 3.1 ± 0.2 GcV~ 



which agrees with the result for D71" (|24|) . This can be 
compared to a calculation using heavy meson ChPT [TtJ 
in which clo/hd-k ~ 10.5 GeV -2 is obtained. Just like 
similar heavy meson ChPT calculations for the Dit scat- 
tering length [111, [lH [ll| this value is somewhat lower 
than our result. 

With our restricted number of phase shift points a 
general extraction of resonance parameters for two res- 
onances is not possible. We therefore have to appeal 
to knowledge from the heavy-quark limit Q , outlined in 
the introduction, in order to extract anything more than 
values of the S-wave phase shift in this channel. One 
expects two J p = 1+ resonances in this limit @: one 
broad resonance with j p = ^ + which decays into D*tt in 

S-wave, and one narrow resonance with j p = | + which 
only decays to D*ir in D-wave. This qualitatively agrees 
with the experiment, where the Di(2430) is broad with 
T = 384±^° MeV, while the Di(2420) is fairly narrow 
with T = 27.1(2.7) MeV 0. The presence of addi- 
tional levels in Fig. [7] is related to resonances, and we 
will assume that the energy level unaffected by the in- 
clusion of meson-meson interpolating fields (red boxes in 
Figurc[5]) corresponds to the narrow Z?i(2420). This res- 
onance is below D-wave threshold D*(1)tt(— 1) for our 
L, so £>i(2420) in our simulation is stable against the 
only decay channel in the heavy quark limit. Away from 
this limit, this state might also develop a small width in 
S-wave, but we expect that the £>i(2420) in our simula- 
tion is even narrower than in experiment. It is therefore 
assumed that this narrow Di(2420) docs not affect the 
D*ir phase shifts at the positions of the other three lev- 
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els, which arc then dominated by the broad resonance 
Z?i(2430). We therefore drop the phase-shift point from 
level n = 2 and proceed to extract a resonance mass and 
width for the £>i(2430) by a linear Breit-Wigncr fit Q 
of the remaining three levels n = 1, 3, 4 in Table |VnT] 
and Figure [9] Under this assumption 



fl ^ = 2.01(15)(02)GeV, 



/7? 



lat 
-Di 



M— 
lvl is 



380.7(20.0) (4.0)MeV, 



(28) 
(29) 



with A% = \{M D + 3M D *). The mass difference with 
respect to the spin-averaged IS state is plotted in Fig- 
ure [lOl The coupling can again be compared to the up- 
per limit from the total width of the Z?i(2430) given by 
9diDtt — 2.50(40). Considering the statistic and system- 
atic uncertainties of our simulation both resonance pa- 
rameters are in reasonable agreement with experiment. 
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FIG. 10. Energy differences AE = E - \ (M D + 3M D *) for 
D meson states on the lattice and in experiment; reference 
spin-averaged mass is \{Md + 3Md*) ~ 1971 MeV in exper- 
iment. Magenta diamonds give resonance masses for states 
treated as resonances in the present simulation. Energy levels 
as extracted in a finite box are given by blue crosses (present 
simulation) and red stars (simulation ||). Established ex- 
perimental states are depicted with black lines or gray boxes 
with a solid black outline: the height indicates experimental 
uncertainty in the resonance mass. 



3. Compilation of D meson results 

Figure [TU] summarizes our results for the low-lying D 
meson states compared to well-established experimental 
results [lH (black lines and boxes). Only the masses (ma- 
genta diamonds) of states corresponding to the broad res- 
onances Z?5(2400) and Z?i(2430) have been extracted tak- 
ing their resonance nature into account explicitly. These 
are the only two states (among six 15 and IP states) 
that are expected to be broad in the m c —¥ oo limit and 
are broad in experiment. The remaining four states can 



Mass difference 


This paper [MeV] 


Experiment [MeV] 


L>5(2400) 


350.8 ±20.2 ±3.7 


347 ± 29 


£>i(2430) 


380.7 ±20.0 ±4.0 


456 ± 40 


£>i(2420) 


448.7 ± 14.1 ±4.7 


449.9 ± 0.6 


D2(2460) 


508.2 ± 17.1 ± 5.3 


491.2 ±0.7 


D(2550) 


599.5 ±24.7 ±6.3 


568.0 ±8.2 


D* (2600) 


674.5 ±28.0 ±7.1 


637.3 ± 3.5 


D(2750) 


793.3 ± 16.4 ±8.3 


781.0 ±3.2 


7-1*/ 

D 


745.8 ± 23.2 ± 7.8 




2S-TS 


655.8 ± 24.6 ± 6.9 


619.9 ±3.7 


IS hyperfine 


129.4 ± 1.8 ± 1.4 


140.65 ±0.10 


2S hyperfine 


75.0 ± 26.9 ± 0.8 


69.3 ± 8.9 



TABLE IX. Mass differences in the D meson spectrum, com- 
pared to experimental values (well established states calcu- 
lated from [33| and others from |(|). Bars denote spin- 
averaged values. For the results of this paper, the first error 
denotes the statistical uncertainty and the second error de- 
notes the uncertainty from setting the lattice scale. Regarding 
the scale-setting uncertainty similar remarks to the Charmo- 
nium case apply. In addition there is a non- negligible error 
from the uncertainty in the determination of k c for all spin- 
dependent quantities. It should be stressed that the gauge 
ensembles at our disposal do not allow for a continuum and 
infinite volume extrapolation. Consequently qualitative but 
not quantitative agreement is expected. 



be treated as stable on our lattice, so we equate their 
masses to the energy levels determined from correlation 
functions using only qq interpolators. We make the same 
assumption for the states in the J p = 2~ channel and for 
excited states with J p = 0~, 1~. For the hadronically 
stable states (D, D* at our simulated pion mass) neglect 
of explicit coupling to multi-hadron states should be a 
good approximation. For narrow states above hadronic 
thresholds one might expect the neglect of explicit cou- 
pling to result in a mass shift comparable to the hadronic 
width. In addition to the results from this work we also 
display the results from [H as red stars. In this case 
masses of all states correspond to energy levels deter- 
mined directly from correlation functions using only qq 
interpolators. As already observed for charmonium, our 
results for the IP and 2S states come out at somewhat 
lower mass than in the previous simulation 0. As we 
are working with a slightly improved heavy-quark treat- 
ment, different sources, a different pion mass, a different 
volume and a different scale setting procedure it is not 
clear what exactly causes this difference. 

In addition to the well established states some new 
resonances were recently observed by the BaBar collab- 
oration Id] . In particular BaBar observes two new res- 
onance^ -D(2550) and £>*(2600) which are interpreted 



For all D meson results we always compare to the neutral states. 
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in the literature as the 2S states [64U69I ]. In addi- 
tion there is also evidence for a state at a mass of 
2752.4 ± 1.7 ± 2.7 MeV in D*ir and an observation of 
a state at a mass of 2763.3 ± 2.3 ± 2.3 MeV in Dir. While 
these signals are interpreted as a single state in the PDG 
[33d and by some authors E3 , others [66T - I691 ] prefer 
the interpretation as two different states. If interpreted 
as two different states, the D(2750) is most commonly in- 
terpreted as either the lowest or first excited state in the 
2" channel, the D 2 or D' 2 . The L»*(2760) is most com- 
monly interpreted as the ground state in the J p = 3~ 
channel, the D^, although quantum numbers J p = 1~ 
can not be excluded either, especially if the two observed 
signals come from the same resonance. 

In Figure [TT1 our results are plotted again with the most 
likely assignments of these new BaBar states under the 
assumption of two distinct resonances for the Z?(2750) 
and .D*(2760) signals. Table HXl provides the numerical 
values for the mass splittings and also includes our re- 
sults for the IS and 2S hypcrfinc splittings and for the 
2S— IS splitting, where correlations have been taken into 
account. 
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FIG. 11. Energy differences AE = E - \{Md + 3M D « ) for D 
meson states in the present simulation and in experiment; the 
reference spin-averaged mass is \{Md + 3Md*) ~ 1971 MeV 
in experiment. Magenta diamonds give resonance masses for 
states treated as resonances in the present simulation. Masses 
extracted as energy levels in a finite box are displayed as blue 
crosses. Established experimental states are depicted with 
black lines or gray boxes with a solid black outline: the height 
indicates experimental uncertainty in the resonance mass. In 
addition to these well-established states the plot also shows 
energy levels from a recent publication by the BaBar collabo- 
ration @ as green boxes with a dotted black outline, choosing 
a set of possible quantum number assignments which seems 
to be favored in the literature [6414691 ]. For further comments 
regarding this assignment please refer to the text. 



V. SUMMARY & DISCUSSION 

Charm-light mesons were studied using a dynamical 
lattice QCD simulation with two flavors of light quarks. 
The main goal was the treatment of the experimentally 
broad scalar and axial D mesons as resonances in Dtt and 
D*ir scattering. The distillation method was used for the 
calculation of quark propagation and this facilitates the 
construction of the correlators that incorporate Dtt and 
D*tt operators in addition to the usual quark-antiquark 
ones. 

The heavy quark was treated using the Fermilab ap- 
proach and the charm quark mass was tuned to fit the 
kinetic mass of the spin-averaged IS* charmonium states. 
As a check, the kinetic masses for spin-averaged S-wave D 
and D s mesons were also calculated. At our final choice 
k c = 0.123 the tuned kinetic masses agree with experi- 
mental values to better that 2%. 

The low-lying charmonium spectrum was calculated 
first, to validate our heavy-quark treatment. The dis- 
tillation method combined with a large basis of quark- 
antiquark operators allowed the extraction of the ground 
and a number of excited states. The resulting spectrum 
for various J PC up to J = 3 in Fig. [5] shows overall 
good agreement with experiment. An interesting feature 
of the calculation is the observation of a state in the 1 ++ 
channel very close in mass to the V(3872). Wc however 
stress that within our systematic uncertainty we can not 
rule out the possibility of quantum numbers 2 ^ for the 
X(3872). 

Encouraged by that, the S-wave phase shifts were cal- 
culated for Dtt scattering with J p = + and D*ir scat- 
tering with J p = 1 + , focusing on isospin 1/2 channels 
where resonances appear. Following the Luscher method, 
we first extracted the discrete energies of the £>W 7r sys- 
tem with total zero momentum in a finite box. The 
energy levels were obtained using quark-antiquark and 
two-meson operators in the correlation functions. The 
Liischcr formula then renders the phase shift for levels in 
the elastic region. 

In the J p = + channel, we extract values of the phase 
shift at three different relative momenta. Only one low- 
lying resonance is expected in this channel and, assuming 
a Brcit-Wigner shape, a resonance mass and width were 
extracted. The resulting resonance mass is 351 ±21 MeV 
above the spin-average j(nijj + 3m£>*). This agrees with 
the mass of the observed resonance -Dp (2400), whose 
mass is 347 ± 29 MeV above the spin average in experi- 
ment. We parametrized the width T = g 2 p* /s, and the 
resulting Dq -> Dn coupling g lat = 2.55 ± 0.21 GcV is 
close to the experimental value g exp < 1.92 ± 0.14 GeV. 

The J p — 1 + channel is more complicated due to the 
presence of two nearby axial resonances. Four energy lev- 
els were extracted. One of the levels was essentially un- 
affected whether the D*ir interpolators were included or 
not, so wc associated this level with the narrow £q(2420). 
The remaining energy levels were used in a Breit-Wigncr 
fit to obtain resonance parameters which are associ- 



14 



ated with the broad Z?i(2430). The resonance mass is 
found at 381 ± 20 MeV above j(mu + 3m£>»), which is 
slightly lower than the experimental value 456 ± 40 MeV, 
while the coupling g lat = 2.01 ± 0.15 GeV agrees with 
g ex P < 2.50 ± 0.40 GeV. 

The main results for the D meson spectrum are com- 
piled in Fig. 1111 where resonance masses for scalar and 
axial mesons arc shown together with our results for other 
J p . The later were calculated using just quark-antiquark 
operators and by equating the masses to the energy lev- 
els. Overall good agreement with experimental values of 
the well established states was obtained. Furthermore 
additional energy levels were observed in the vicinity of 
some of the resonances discovered recently by the BaBar 
collaboration jf|. 

In addition to the resonance parameters, the S-wave 
scattering lengths do were determined from the ground 

1—1/2 

states. The resulting a = 0.81 ± 0.14 fm for Dir and 

ag =1/2 = 0.81 ± 0.17 fm for D*n apply for = 266 
MeV in our simulation. 

The experimental observation of new D and D s meson 
states in the past decade led to a number of challenges 
for theory and new ideas emerged. An example is the 
suggestion that explicit ss content should be included in 
the structure of some charm-light meson states (see for 
example [70| and references therein). In this context, it 
is interesting that the present simulation results in favor- 
able agreement with experiment without the inclusion of 
strange quark content. 



Appendix A: Interpolating fields 

1. Charmonium 

Table |X] lists the interpolating fields used in Section 
IIIII to study the low-lying charmonium spectrum. The 
symbol V; is used for a single covariant derivative. For a 
Laplacian-like structure we use A = X^=i V,Vi. In ad- 
dition we also consider structures with a symmetrized 
second derivative of type £>; = \euk jVj 'Vfc first pro- 
posed in [7lj and previously used in [72j . Summation 
over repeated roman indices is implied and e^fc denotes 
the Levi-Civita symbol in three dimensions. For inter- 
polating fields in the E representation some non-trivial 
Clebsch-Gordan coefficients are needed. They are 





Qm — 


1 

71' 


Q122 = - 


1 

71' 


1 

71' 


Q222 = 


1 

"71 


Q233 = 


tV (A1) 


i 

71' 


Q122 = 


1 

"71 


Q1.33 = 


2 

71' 




Q211 


1 

"71 


Q222 = 


1 

71' 



Qm 



For completeness, Table IXII lists the interpolator 
choice, timeslice ioj fit range and type as well as the fit re- 
sults and the x 2 /d.o.f for all fits performed to determine 
energy levels related to charmonium states. 
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TABLE X: Table of cc interpolators used for charmonium 
states. They are sorted by irreducible representation of 
the octahedral group Oh and by quantum numbers PC. 
The reduced lattice symmetry implies an infinite number 
of continuum spins in each irreducible representation of 
the octahedral group. For operators, repeated roman in- 
dices indicate summation. The quantity "ft denotes the 
Dirac matrix for the time direction. 



Lattice 


Quantum numbers 


Interpolator 


Operator 


irrep 


J FC in irrep 


label 




A1-+ 


0-+ 4"+,... 


1 


97517 






2 


qjt-ysq 






3 


97*7^75 V 'tq 






4 


g^7sVig 






5 


g^i7t7sVtiJ 






6 


qA'ysAq 






7 


<r 

gA7t75 Aq 






8 


<r- — > 

qAjt-yiJsViq 


A1++ 


++ ,4++, ... 


1 


qq 






2 


qiNiq 






3 


qjt'fiViq 






4 





Continued in next column 
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TABLE X — continued from previous column 



Lattice 


Quantum numbers 


Interpolator 


Operator 


irrep 


T pfi . 

J m irrep 


label 
5 
6 
7 


g"AAg 
gA7iVig 
gA7t7iVig 


T~ 


I", 3—, 4", ... 


1 

2 
3 
4 
5 
6 
7 
8 
9 
1U 
11 
12 
13 
14 


<77;<7 
<77t7i<? 
gV^g 

q^~i^fiV^q 
qVi^tlNiq 

gA7i Ag 

tr ^ 
q/X-ytH Ag 

gA~vtg 

qAeij k jjj 5 S7 k q 

q\ajkhj D kq 

q\eijkhtJjD k q 
qibBiq 


T++ 


1++, 3++, 4++, ... 


1 

2 
3 

4 

5 
6 
7 
8 


qiUbq 
qtijkij^kq 

gey fe 7i7jV/bg 
g^i7i75 V^g 
g A 7^75 A g 
gAe i jfc7jVj I g 
qAe ij kjt'yjVkq 

q\ti ] k\'ys'y]D k q 


T+~ 


l + ~ 3 + ~ 4+~ ... 


1 

2 
3 

4 

5 

D 

7 
8 


g7t7i75<7 
975 Vig 
g7t75Vig 
5^i7t7i75Vig 
gA7t7;75 Ag 
g A75 V^g 
gA7i7 5 Vig 
g|eijfe|7t7573- D fc9 




2 ,3 ,4 


1 

O 
z 


9 s |eyfc|7j75?D fc 
9^l e »jfc|7j759i5 fc 


T -+ 

± 2 


2" + , 3"+, 4"+, ... 


1 

Z 

3 
4 


g|eijfe|7t7i75V fc g 
9 A |ejjfe |7t7j75 Vfeg 
975 A? 


rp+ + 


2++, 3++, 4++, ... 


1 

2 


q \Ujk\7jVkq 
q tijfc 7*7jVfcg 


Continued in next column 



TABLE X — continued from previous column 



Lattice 


Quantum numbers 


Interpolator 


Operator 


irrep 


T pfi . 

(/ in niep 


IdUcl 








Q 
O 


/-r /\ . . r r\i r \ / r ft 

qL\ \ ttjk \j3 v fcy 






4 








K 





— 


9 4 — 

Z ,41 , . . . 


\ 


fit ) • 'I *V ■ ^Vr- \ / 1 / J 

yVajfc7j75 v fcy 






9 
z 


<? A Q y fc 7j 75 V fc g 






Q 
O 


q^iijklj^kq 






4 


H^Cijk ft IJ^kH 




0-+ 4-+ 

z , t , . . . 


1 








9 
z 


gvijfc7*7j75 v fcy 




2 ,4 , . . . 


1 

1 


qQijkJj Vfcg 






9 
z 


qQijk'Yt'Yj'V kq 






■3 
O 


g^^rzjfc7j v fcy 






4 


H^^iijk ft }j vfcg 






5 


qQ'ijk'ysyjDtq 


At + 


3++, 6++, ... 


1 


qi5liDiq 




3+", 6+",... 


1 


qitisiiDiq 




3—, 6—, ... 


1 


qjiDtq 






2 


qjtJiDiq 



2. Fitting methodology 

Depending on the channel the full interpolator basis 
is pruned to a less noisy subset. We cither fit with a 
single exponential at large time separations or with two 
exponentials starting at smaller time separations. A jack- 
knife estimate of the covariance matrix on the ensem- 
ble average is used to perform correlated fits. To build 
the pseudo-inverse of the matrix we perform a singular 
value decomposition and exclude very tiny singular val- 
ues when the ratio of largest over smallest singular values 
gets close to machine precision. In these cases it is nec- 
essary to remove the corresponding number of degrees of 
freedom from the fit. Table IXII lists our choices for the 
interpolator sets, timeslice to and fit ranges as well as 
the fit results and x 2 /d.o.f. for the fits. The basis used 
is indicated by the interpolator numbers from Table IXl 



3. D mesons 

Tabic IXIll lists the quark- antiquark interpolating fields 
used in Section IIVI The notation is the same than for 
charmonium in the previous section. In addition Table 
IXIIII lists the interpolator choice, timeslice to, fit range 
and type as well as the fit results and the x 2 /d.o.f for 
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channel 


state 


interpolators 


to 


fit range 


fit type 


E n a 


x7d-o.f 


Al~ + 


I 


1 3 4 6 S 

1,0,1,0,0 


2 


3-27 


2 exp 




13.68/7 




2 


1 3 4 6 8 

1,0,1,0,0 


2 


3-18 


2 6Xp 


1 8694(32) 


15.47/7 




■) 
o 


13 4 6 8 
1,0,1,0,0 


9 
— 


3 1 1 
1 1 


2 6Xp 




1 6Q /5 


A1++ 


1 


12 3 4 6 7 

1,Z:,0,1,U, ( 


3 


4-21 


2 6Xp 


1 7475(29) 


4.57/8 




2 


1 2 3 4 6 7 
1,^,0,1,0, 1 


3 


7-12 


1 GXp 


2 021 (1 3) 


0.63/4 




3 


1 2 3 4 6 7 

i,z,, 0,1,0, ( 


3 


7-11 


1 GXp 


2 1 54(31 ) 

ij . 1Ut;( 'Jl / 


1.10/3 


Tl~~ 


1 


1 2 3 4 5 6 7 8 11 12 


2 


3-27 


2 GXp 


1 54171(43) 


6.75/6 




2 


1 9 3 4 5 6 7 8 11 19 


2 


3-16 




1 9058(33) 


18.35 /7 




3 


1 9 3 4 5 6 7 8 11 19 


2 


3-15 


2 exp 


1 9838(47) 


3.80 /7 




4 


1 234 5 67811 12 


2 


3-14 


2 GXp 


2 0144(45) 


7.45/8 




5 


1 234567811 19 


2 


3-9 




9 1 79(1 7 s ; 


2.28/3 




g 


1 234 5 67811 12 


2 


3-9 


2 GXp 


2 189(40) 


0.9/3 


T1++ 


1 


1 3 4 5 7 8 
1,0,1,0, 1 ,0 


3 


4-14 


2 GXp 


1 7921 (26) 


8.27/7 




2 


13 4 5 7 8 
1,0,1,0, 1 ,0 


3 


8-13 


1 GXp 


2 041(18') 


0.64/4 




3 


1 3 4 5 7 8 
1,0,1,0, ( ,0 


3 


7-12 




9 1 46(1 5) 


2.06/4 


Tl + ~ 


1 


1 3 4 5 7 8 

1,0,1,0, / ,0 


3 


4-14 


2 GXp 


1 7998(32) 


6.14/7 




2 


1 3 4 5 7 8 

1,0,1,0, / ,0 


3 


8-12 


1 GXp 


2 053(20) 


1.57/3 




3 


1 3 4 5 7 8 

1,0,1,0, ( ,0 


3 


7-11 




9 1 45(1 7) 


1.21/3 


T2~~ 


1 


1 2 
1,- 


2 


3-16 


2 GXp 


2 0031 (43) 


4.09/8 




2 


1 9 
1.^ 


2 


3-14 


2 GXp 


2 063(1 3) 


5.37/7 


T2~ + 


1 


12 3 4 
1,^,0,1 


2 


3-14 


2 GXp 


2 0001 (41 ) 


8.05/7 




2 


19 3 4 
i,z,o,i 


2 


7-10 




9 190(56) 


1.01/2 


T2 ++ 


1 


12 3 5 


2 


4-17 


2 GXp 


1 81 90(31 ) 


5.48/7 




2 


19 3 5 


2 


8-13 




9 061 (95) 


0.76/4 




3 


1,2,3,5 


2 


7-12 


1 GXp 


2.143(11) 


1.53/4 


A2— 


1 


1,2 


2 


3-14 


2 GXp 


2.0124(38) 


2.61/7 


A2++ 


1 


1 


2 


7-12 


1 GXp 


2.164(16) 


1.85/4 


A2+- 


1 


1 


2 


7-12 


1 GXp 


2.155(20) 


2.04./4 



TABLE XI. Fit details for the charmonium states. Interpolators are listed in Table |X] 



those D meson-levels which are either stable under the strong interaction or are narrow and have not been prop- 
erly treated as resonances in our current study. 
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Lattice 


Quantum numbers 


Interpolator 


Operator 


irrep 


J m irrep 


label 




AY~ 


0" 4" 


1 


- / 






2 








3 


!77t7i75_Vi<?' 






4 


57*76 Vi^' 






5 








6 


g^7t7sVtg' 


A1+ 


0+,4+, ... 


1 








2 








3 








4 


g^vtg' 


Tr 


1- 3-, 4-, ... 


1 


974 ?' 






2 


*i / 1 (If 






3 


gvtg' 






4 


q^ijkJjJs'V kq' 






5 


g7*V;g' 






6 


qeijfe7«7i75Vfcg' 






7 


gyi7iVig' 






8 


^7t7iVt'7' 


T+ 


1+ 3+ 4+ ... 


1 








2 


y <■ j re /J w ^ 






3 


oe« .-{■'Yi'YiVi-o' 






4 


97t7i75?' 






5 


qysVjtf 






6 


?7t75Vig' 






7 


g^7i75vtg' 






8 


g^7i7i75Vig' 




2", 3", 4", ... 


1 


g s eijfc 7j7sVfeg' 






2 


g|eijfe|7t7i75V fc g' 


n 


2+,3+,4+, ... 


1 


g|e»jfc 7jVfcg' 






2 


g|e lJ fe|7i7jV fe g' 



TABLE XII. Table of interpolators used for D mesons; in 
addition we use Di (|21|) and D*7r ([23} interpolators for irreps 
and T+. Interpolators are sorted by irreducible represen- 
tation of the octahedral group Oh and by the parity quantum 
number P. The reduced lattice symmetry implies an infinite 
number of continuum spins in each irreducible representation 
of Oh- For operators, repeated roman indices indicate sum- 
mation. The quantity jt denotes the Dirac matrix for the 
time direction. 
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channel 


state 


interpolators 


to 


fit range 


fit type 


E n a 


X 2 /d.o.f 


Al- 


1 


1,2,5 


2 


4-28 


2 exp 


0.97952(78) 


15.63/14 




2 


1,2,5 


2 


8-14 


1 exp 


1.417(16) 


0.68/5 


Tl~ 


1 


1,2,3,7 


2 


4-28 


2 exp 


1.0608(15) 


12.45/14 
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1,2,3,7 
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8-12 


1 exp 


1.464(18) 


2.63/3 
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1,2,3,7 
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8-12 


1 exp 


1.509(15) 


0.56/3 


T2~ 
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1,2 
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7-11 


1 exp 


1.539(10) 


0.40/3 




2 


1,2 


5 


7-11 


1 exp 


1.557(15) 


3.10/3 


T2+ 


1 


1,2 


3 


4-15 


2 exp 


1.360(11) 


4.85/8 



TABLE XIII. Fit details for the D-meson channels that incorporate only uc interpolators, listed in Table IXlTl 
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